function pt
global m d e
m=2;b=[0,0.3,0.2];p=[0,0,1];
px=[4.6;4.5;4.5];py=[3.5;1.8;0.4];
strdd{1}='无阻尼';
strdd{2}='阻力正比于v';
strdd{3}='阻力正比于v^2';
figure
for i=1:3
    d=b(i); e=p(i);
    [t,y]=ode45(@ptfun,[0:0.01:10],[0,3,0,5]);
    H{i}=max(y(:,3))
    T{i}=t(find(y(:,3)==H{i}))
    vx0{i}=y(end,2)
    subplot(2,1,1)
    axis([0 6 -70 2]);
    hold on
    xlabel('x');
    ylabel('y');
    comet(y(:,1),y(:,3));
    subplot(2,1,2)
    axis([0 10 0 4])
    hold on
    xlabel('t')
    ylabel('dx/dt')
    text(px(i),py(i),strdd{i});
    comet(t,y(:,2))
end
function ydot=ptfun(t,y)
    global m d e
    ydot=[y(2);-d/m*y(2)*(y(2).^2+y(4).^2)^(e/2);y(4);
        -9.8-d/m*y(4)*(y(2).^2+y(4).^2)^(e/2)];